Importing the Data¶
In [1]:
# @title
import pandas as pd
from sklearn.model_selection import cross_val_score, TimeSeriesSplit
data = pd.read_csv(
'https://raw.githubusercontent.com/Garve/datasets/4576d323bf2b66c906d5130d686245ad205505cf/mmm.csv',
parse_dates=['Date'],
index_col='Date'
)
X = data.drop(columns=['Sales'])
y = data['Sales']
Parameter Values¶
In [2]:
# @title
N = 200 # number of observations
# previous marketing mix modeling has given us these parameters
tv_coef = 10000 # α
tv_lags = 4 # ℓ
tv_carryover = 0.5 # λ
tv_saturation = 0.002 # β
radio_coef = 8000
radio_lags = 2
radio_carryover = 0.2
radio_saturation = 0.0001
banners_coef = 14000
banners_lags = 0
banners_carryover = 0.2
banners_saturation = 0.001
Carryover Matrices¶
In [3]:
# @title
from scipy.linalg import toeplitz
import numpy as np
# Much cleaner!
tv_carryover_matrix = toeplitz(
tv_carryover**np.arange(N), # first column
np.zeros(N) # first row (all zeros = lower triangular)
)
radio_carryover_matrix = toeplitz(
radio_carryover**np.arange(N),
np.zeros(N)
)
banners_carryover_matrix = np.eye(N)
Prediction¶
In [4]:
# @title
# Calculate predicted sales per period (corrected formula)
sales_per_period = (
tv_coef * (1 - np.exp(-tv_saturation * tv_carryover_matrix @ data["TV"]))
+ radio_coef * (1 - np.exp(-radio_saturation * radio_carryover_matrix @ data["Radio"]))
+ banners_coef * (1 - np.exp(-banners_saturation * banners_carryover_matrix @ data["Banners"]))
)
# Total sales comparison
total_actual_sales = y.sum()
total_predicted_sales = sales_per_period.sum()
print(f"Total Actual Sales: ${total_actual_sales:,.2f}")
print(f"Total Predicted Sales: ${total_predicted_sales:,.2f}")
print(f"Difference: ${total_predicted_sales - total_actual_sales:,.2f}")
print(f"Percentage Difference: {((total_predicted_sales - total_actual_sales) / total_actual_sales * 100):.2f}%")
Total Actual Sales: $2,133,628.30 Total Predicted Sales: $3,847,201.38 Difference: $1,713,573.08 Percentage Difference: 80.31%
Prediction using Optimized Parameters¶
In [5]:
# @title
from scipy.optimize import minimize
from scipy.linalg import toeplitz
import numpy as np
# Prepare data
X_tv = data["TV"].values
X_radio = data["Radio"].values
X_banners = data["Banners"].values
y_actual = y.values
N = len(y)
def mmm_prediction(params, X_tv, X_radio, X_banners, N):
"""Generate sales prediction given parameters"""
tv_coef, tv_carryover, tv_saturation = params[0:3]
radio_coef, radio_carryover, radio_saturation = params[3:6]
banners_coef, banners_saturation = params[6:8]
# Create carryover matrices with fixed lags
tv_lags, radio_lags = 4, 2
tv_matrix = toeplitz(tv_carryover**np.arange(N), np.zeros(N))
radio_matrix = toeplitz(radio_carryover**np.arange(N), np.zeros(N))
banners_matrix = np.eye(N)
# Calculate contributions
tv_contribution = tv_coef * (1 - np.exp(-tv_saturation * tv_matrix @ X_tv))
radio_contribution = radio_coef * (1 - np.exp(-radio_saturation * radio_matrix @ X_radio))
banners_contribution = banners_coef * (1 - np.exp(-banners_saturation * banners_matrix @ X_banners))
return tv_contribution + radio_contribution + banners_contribution
def objective(params, X_tv, X_radio, X_banners, y_actual, N):
"""Mean Absolute Percentage Error"""
y_pred = mmm_prediction(params, X_tv, X_radio, X_banners, N)
mape = np.mean(np.abs((y_actual - y_pred) / y_actual)) * 100
return mape
# Initial guess (your original parameters)
initial_params = [
10000, 0.5, 0.002, # TV: coef, carryover, saturation
8000, 0.2, 0.0001, # Radio: coef, carryover, saturation
14000, 0.001 # Banners: coef, saturation
]
# Bounds to keep parameters reasonable
bounds = [
(0, 50000), # tv_coef
(0, 0.99), # tv_carryover
(0.0001, 0.01), # tv_saturation
(0, 50000), # radio_coef
(0, 0.99), # radio_carryover
(0.0001, 0.01), # radio_saturation
(0, 50000), # banners_coef
(0.0001, 0.01), # banners_saturation
]
print("Optimizing parameters...")
result = minimize(
objective,
initial_params,
args=(X_tv, X_radio, X_banners, y_actual, N),
method='L-BFGS-B',
bounds=bounds,
options={'maxiter': 1000}
)
# Extract optimized parameters
opt_params = result.x
print(f"\nOptimization successful: {result.success}")
print(f"Final MAPE: {result.fun:.2f}%")
print("\nOptimized Parameters:")
print(f"TV: coef={opt_params[0]:.2f}, carryover={opt_params[1]:.4f}, saturation={opt_params[2]:.6f}")
print(f"Radio: coef={opt_params[3]:.2f}, carryover={opt_params[4]:.4f}, saturation={opt_params[5]:.6f}")
print(f"Banners: coef={opt_params[6]:.2f}, saturation={opt_params[7]:.6f}")
# Compare results
y_pred_optimized = mmm_prediction(opt_params, X_tv, X_radio, X_banners, N)
print("\n=== Total Sales Comparison ===")
print(f"Actual: ${y_actual.sum():,.2f}")
print(f"Original Model: ${sales_per_period.sum():,.2f}")
print(f"Optimized Model: ${y_pred_optimized.sum():,.2f}")
Optimizing parameters... Optimization successful: True Final MAPE: 12.19% Optimized Parameters: TV: coef=10000.00, carryover=0.7257, saturation=0.000100 Radio: coef=8000.00, carryover=0.2909, saturation=0.000104 Banners: coef=14000.00, saturation=0.000145 === Total Sales Comparison === Actual: $2,133,628.30 Original Model: $3,847,201.38 Optimized Model: $2,085,086.69
In [ ]:
# @title
import matplotlib.pyplot as plt
fig, axes = plt.subplots(2, 1, figsize=(14, 8))
# Time series comparison
axes[0].plot(data.index, y_actual, label='Actual', linewidth=2)
axes[0].plot(data.index, y_pred_optimized, label='Predicted', linewidth=2, alpha=0.7)
axes[0].set_title('Actual vs Predicted Sales Over Time')
axes[0].legend()
axes[0].set_ylabel('Sales')
# Residuals
residuals = y_actual - y_pred_optimized
axes[1].plot(data.index, residuals)
axes[1].axhline(y=0, color='r', linestyle='--')
axes[1].set_title('Residuals')
axes[1].set_ylabel('Residual')
plt.tight_layout()
plt.show()
Contributions by channel¶
In [ ]:
# @title
# Recreate the carryover matrices with optimized parameters
tv_matrix = toeplitz(opt_params[1]**np.arange(N), np.zeros(N))
radio_matrix = toeplitz(opt_params[4]**np.arange(N), np.zeros(N))
banners_matrix = np.eye(N)
# Calculate individual channel contributions
tv_contribution = opt_params[0] * (1 - np.exp(-opt_params[2] * tv_matrix @ X_tv))
radio_contribution = opt_params[3] * (1 - np.exp(-opt_params[5] * radio_matrix @ X_radio))
banners_contribution = opt_params[6] * (1 - np.exp(-opt_params[7] * banners_matrix @ X_banners))
print("Total Sales Contribution by Channel:")
print(f"TV: ${tv_contribution.sum():,.2f} ({tv_contribution.sum()/y_pred_optimized.sum()*100:.1f}%)")
print(f"Radio: ${radio_contribution.sum():,.2f} ({radio_contribution.sum()/y_pred_optimized.sum()*100:.1f}%)")
print(f"Banners: ${banners_contribution.sum():,.2f} ({banners_contribution.sum()/y_pred_optimized.sum()*100:.1f}%)")
Total Sales Contribution by Channel: TV: $1,146,349.83 (55.0%) Radio: $403,195.47 (19.3%) Banners: $535,541.40 (25.7%)
Goodness of Fit¶
In [ ]:
# @title
# Calculate R-squared
ss_res = np.sum((y_actual - y_pred_optimized)**2) # Residual sum of squares
ss_tot = np.sum((y_actual - y_actual.mean())**2) # Total sum of squares
r_squared = 1 - (ss_res / ss_tot)
print("\n=== Model Performance Metrics ===")
print(f"R-squared: {r_squared:.4f}")
print(f"Adjusted R-squared: {1 - (1-r_squared)*(N-1)/(N-8-1):.4f}") # 8 parameters
print(f"MAPE: {result.fun:.2f}%")
print(f"RMSE: ${np.sqrt(ss_res/N):,.2f}")
# Alternative: using sklearn
from sklearn.metrics import r2_score, mean_absolute_percentage_error, mean_squared_error
print("\n=== Using sklearn ===")
print(f"R²: {r2_score(y_actual, y_pred_optimized):.4f}")
print(f"MAPE: {mean_absolute_percentage_error(y_actual, y_pred_optimized)*100:.2f}%")
print(f"RMSE: ${np.sqrt(mean_squared_error(y_actual, y_pred_optimized)):,.2f}")
=== Model Performance Metrics === R-squared: 0.7101 Adjusted R-squared: 0.6979 MAPE: 12.19% RMSE: $1,450.54 === Using sklearn === R²: 0.7101 MAPE: 12.19% RMSE: $1,450.54
Optimal Budget Allocation Using CVXPY¶
In [ ]:
# @title
import cvxpy as cp
import numpy as np
from scipy.linalg import toeplitz
# Extract optimized parameters
tv_coef = opt_params[0]
tv_carryover = opt_params[1]
tv_saturation = opt_params[2]
radio_coef = opt_params[3]
radio_carryover = opt_params[4]
radio_saturation = opt_params[5]
banners_coef = opt_params[6]
banners_saturation = opt_params[7]
# Create carryover matrices with optimized parameters
tv_carryover_matrix = toeplitz(tv_carryover**np.arange(N), np.zeros(N))
radio_carryover_matrix = toeplitz(radio_carryover**np.arange(N), np.zeros(N))
banners_carryover_matrix = np.eye(N)
# Get original total budget
original_total_spends = data[["TV", "Radio", "Banners"]].sum().sum()
# Declare CVXPY variables (N=200 per channel)
tv = cp.Variable(N)
radio = cp.Variable(N)
banners = cp.Variable(N)
# Constraints
constraints = [
tv >= 0,
radio >= 0,
banners >= 0,
cp.sum(tv + radio + banners) <= original_total_spends,
]
# Objective function using optimized parameters
objective = cp.Maximize(
tv_coef * cp.sum(1 - cp.exp(-tv_saturation * tv_carryover_matrix @ tv))
+ radio_coef * cp.sum(1 - cp.exp(-radio_saturation * radio_carryover_matrix @ radio))
+ banners_coef * cp.sum(1 - cp.exp(-banners_saturation * banners_carryover_matrix @ banners))
)
# Create and solve problem (CVXPY auto-selects solver)
problem = cp.Problem(objective, constraints)
problem.solve(verbose=True)
# Display results
print("\n" + "="*60)
print("OPTIMIZATION RESULTS")
print("="*60)
print(f"Solver used: {problem.solver_stats.solver_name}")
print(f"Status: {problem.status}")
print(f"\nOptimal Total Sales: ${problem.value:,.2f}")
print(f"\n--- Budget Allocation ---")
print(f"Original Total Budget: ${original_total_spends:,.2f}")
print(f"Optimized Total Budget: ${(tv.value.sum() + radio.value.sum() + banners.value.sum()):,.2f}")
print(f"\n--- Spend per Channel ---")
print(f"TV: ${tv.value.sum():,.2f} (Original: ${data['TV'].sum():,.2f})")
print(f"Radio: ${radio.value.sum():,.2f} (Original: ${data['Radio'].sum():,.2f})")
print(f"Banners: ${banners.value.sum():,.2f} (Original: ${data['Banners'].sum():,.2f})")
print(f"\n--- Sales Improvement ---")
print(f"Actual Sales (from data): ${y_actual.sum():,.2f}")
print(f"Optimized Sales: ${problem.value:,.2f}")
print(f"Improvement: ${problem.value - y_actual.sum():,.2f}")
print(f"Improvement %: {((problem.value - y_actual.sum()) / y_actual.sum() * 100):.2f}%")
(CVXPY) Jan 17 01:28:50 PM: Your problem has 600 variables, 601 constraints, and 0 parameters. (CVXPY) Jan 17 01:28:50 PM: It is compliant with the following grammars: DCP, DQCP (CVXPY) Jan 17 01:28:50 PM: (If you need to solve this problem multiple times, but with different data, consider using parameters.) (CVXPY) Jan 17 01:28:50 PM: CVXPY will first compile your problem; then, it will invoke a numerical solver to obtain a solution. (CVXPY) Jan 17 01:28:50 PM: Your problem is compiled with the CPP canonicalization backend. (CVXPY) Jan 17 01:28:50 PM: Compiling problem (target solver=CLARABEL). (CVXPY) Jan 17 01:28:50 PM: Reduction chain: FlipObjective -> Dcp2Cone -> CvxAttr2Constr -> ConeMatrixStuffing -> CLARABEL (CVXPY) Jan 17 01:28:50 PM: Applying reduction FlipObjective (CVXPY) Jan 17 01:28:50 PM: Applying reduction Dcp2Cone (CVXPY) Jan 17 01:28:50 PM: Applying reduction CvxAttr2Constr (CVXPY) Jan 17 01:28:50 PM: Applying reduction ConeMatrixStuffing (CVXPY) Jan 17 01:28:50 PM: Applying reduction CLARABEL (CVXPY) Jan 17 01:28:50 PM: Finished problem compilation (took 1.177e-01 seconds). (CVXPY) Jan 17 01:28:50 PM: Invoking solver CLARABEL to obtain a solution.
===============================================================================
CVXPY
v1.6.7
===============================================================================
-------------------------------------------------------------------------------
Compilation
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
Numerical solver
-------------------------------------------------------------------------------
-------------------------------------------------------------
Clarabel.rs v0.11.1 - Clever Acronym
(c) Paul Goulart
University of Oxford, 2022
-------------------------------------------------------------
problem:
variables = 1200
constraints = 2401
nnz(P) = 0
nnz(A) = 42200
cones (total) = 601
: Nonnegative = 1, numel = 601
: Exponential = 600, numel = (3,3,3,3,...,3)
settings:
linear algebra: direct / faer, precision: 64 bit (2 threads)
max iter = 200, time limit = Inf, max step = 0.990
tol_feas = 1.0e-8, tol_gap_abs = 1.0e-8, tol_gap_rel = 1.0e-8,
static reg : on, ϵ1 = 1.0e-8, ϵ2 = 4.9e-32
dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-7
iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
max iter = 10, stop ratio = 5.0
equilibrate: on, min_scale = 1.0e-4, max_scale = 1.0e4
max iter = 10
iter pcost dcost gap pres dres k/t μ step
---------------------------------------------------------------------------------------------
0 +0.0000e+00 -2.3105e+06 2.31e+06 1.00e+00 1.34e+00 1.00e+00 1.00e+00 ------
1 +1.2308e+04 -2.2963e+06 1.88e+02 9.99e-01 1.04e+00 4.08e+00 2.45e-01 7.92e-01
2 +3.0658e+05 -1.9551e+06 7.38e+00 9.77e-01 1.26e+00 7.09e+01 1.38e-02 9.90e-01
3 +1.1475e+06 -9.5952e+05 2.20e+00 9.06e-01 1.09e+00 2.67e+02 3.29e-03 7.71e-01
4 +2.4606e+06 +7.6837e+05 2.20e+00 7.17e-01 6.44e-01 5.87e+02 1.03e-03 7.12e-01
5 +3.4285e+06 +2.6422e+06 2.98e-01 3.23e-01 1.83e-01 6.33e+02 2.44e-04 7.99e-01
6 +3.7013e+06 +3.3954e+06 9.01e-02 1.23e-01 5.85e-02 2.51e+02 7.63e-05 7.92e-01
7 +3.8622e+06 +3.7795e+06 2.19e-02 3.28e-02 1.46e-02 6.98e+01 1.88e-05 7.73e-01
8 +3.9037e+06 +3.8820e+06 5.58e-03 8.55e-03 3.72e-03 1.83e+01 4.80e-06 7.64e-01
9 +3.9143e+06 +3.9094e+06 1.27e-03 1.96e-03 8.48e-04 4.22e+00 1.09e-06 7.92e-01
10 +3.9165e+06 +3.9154e+06 2.82e-04 4.34e-04 1.88e-04 9.36e-01 2.42e-07 7.92e-01
11 +3.9170e+06 +3.9167e+06 6.12e-05 9.43e-05 4.07e-05 2.04e-01 5.25e-08 7.92e-01
12 +3.9170e+06 +3.9170e+06 1.32e-05 2.03e-05 8.78e-06 4.39e-02 1.13e-08 7.92e-01
13 +3.9171e+06 +3.9170e+06 2.85e-06 4.38e-06 1.89e-06 9.47e-03 2.44e-09 7.91e-01
14 +3.9171e+06 +3.9171e+06 6.30e-07 9.71e-07 4.20e-07 2.10e-03 5.41e-10 7.85e-01
15 +3.9171e+06 +3.9171e+06 3.78e-07 5.82e-07 2.52e-07 1.26e-03 3.24e-10 5.07e-01
16 +3.9171e+06 +3.9171e+06 3.51e-07 5.41e-07 2.34e-07 1.17e-03 3.01e-10 2.08e-01
(CVXPY) Jan 17 01:28:51 PM: Problem status: optimal (CVXPY) Jan 17 01:28:51 PM: Optimal value: 2.483e+06 (CVXPY) Jan 17 01:28:51 PM: Compilation took 1.177e-01 seconds (CVXPY) Jan 17 01:28:51 PM: Solver (including time spent in interface) took 6.037e-01 seconds
17 +3.9171e+06 +3.9171e+06 3.38e-07 5.21e-07 2.25e-07 1.13e-03 2.90e-10 1.06e-01
18 +3.9171e+06 +3.9171e+06 3.38e-07 5.21e-07 2.25e-07 1.13e-03 2.90e-10 0.00e+00
19 +3.9171e+06 +3.9171e+06 1.33e-07 2.06e-07 8.88e-08 4.42e-04 1.16e-10 6.34e-01
20 +3.9171e+06 +3.9171e+06 2.86e-08 4.41e-08 1.90e-08 9.46e-05 2.49e-11 7.92e-01
21 +3.9171e+06 +3.9171e+06 1.13e-08 1.75e-08 7.55e-09 3.76e-05 9.85e-12 6.34e-01
22 +3.9171e+06 +3.9171e+06 2.43e-09 3.74e-09 1.62e-09 8.06e-06 2.12e-12 7.92e-01
---------------------------------------------------------------------------------------------
Terminated with status = Solved
solve time = 582.139145ms
-------------------------------------------------------------------------------
Summary
-------------------------------------------------------------------------------
============================================================
OPTIMIZATION RESULTS
============================================================
Solver used: CLARABEL
Status: optimal
Optimal Total Sales: $2,482,938.10
--- Budget Allocation ---
Original Total Budget: $1,336,103.05
Optimized Total Budget: $1,336,103.04
--- Spend per Channel ---
TV: $609,141.90 (Original: $589,241.53)
Radio: $0.08 (Original: $442,717.01)
Banners: $726,961.06 (Original: $304,144.51)
--- Sales Improvement ---
Actual Sales (from data): $2,133,628.30
Optimized Sales: $2,482,938.10
Improvement: $349,309.80
Improvement %: 16.37%
Budget Allocation Optimization with Scipy Minimize (SLSQP)¶
In [ ]:
# @title
from scipy.optimize import minimize
from functools import partial
import numpy as np
from scipy.linalg import toeplitz
# Extract optimized parameters
tv_coef = opt_params[0]
tv_carryover = opt_params[1]
tv_saturation = opt_params[2]
radio_coef = opt_params[3]
radio_carryover = opt_params[4]
radio_saturation = opt_params[5]
banners_coef = opt_params[6]
banners_saturation = opt_params[7]
# Create carryover matrices
N = 200
tv_carryover_matrix = toeplitz(tv_carryover**np.arange(N), np.zeros(N))
radio_carryover_matrix = toeplitz(radio_carryover**np.arange(N), np.zeros(N))
banners_carryover_matrix = np.eye(N)
# Get original budget and channel totals
original_total_budget = data[["TV", "Radio", "Banners"]].sum().sum()
original_channel_spends = data[["TV", "Radio", "Banners"]].sum().values
# Initial guess - original total spends per channel
initial_channel_totals = original_channel_spends
# Bounds - all channel spends must be >= 0
bounds = [(0, None), (0, None), (0, None)]
# Objective function to minimize (we'll negate sales since we want to maximize)
def mmm_objective_by_channel(channel_totals, tv_coef, tv_sat, tv_matrix,
radio_coef, radio_sat, radio_matrix,
banners_coef, banners_sat, banners_matrix,
original_spends, N):
"""
Calculate negative total sales for optimization.
channel_totals: [tv_total, radio_total, banners_total]
Distributes each channel's total proportionally to original daily spends
"""
tv_total, radio_total, banners_total = channel_totals
# Distribute budget proportionally to original daily pattern
tv_spend = original_spends["TV"].values * (tv_total / original_spends["TV"].sum())
radio_spend = original_spends["Radio"].values * (radio_total / original_spends["Radio"].sum())
banners_spend = original_spends["Banners"].values * (banners_total / original_spends["Banners"].sum())
# Calculate contributions with adstock and saturation
tv_contribution = tv_coef * np.sum(1 - np.exp(-tv_sat * tv_matrix @ tv_spend))
radio_contribution = radio_coef * np.sum(1 - np.exp(-radio_sat * radio_matrix @ radio_spend))
banners_contribution = banners_coef * np.sum(1 - np.exp(-banners_sat * banners_matrix @ banners_spend))
total_sales = tv_contribution + radio_contribution + banners_contribution
# Return negative because we're minimizing
return -total_sales
# Budget constraint (equality)
def budget_constraint(channel_totals, budget):
"""Equality constraint: total spend must equal budget"""
return np.sum(channel_totals) - budget
# Create partial function with fixed parameters
partial_objective = partial(
mmm_objective_by_channel,
tv_coef=tv_coef,
tv_sat=tv_saturation,
tv_matrix=tv_carryover_matrix,
radio_coef=radio_coef,
radio_sat=radio_saturation,
radio_matrix=radio_carryover_matrix,
banners_coef=banners_coef,
banners_sat=banners_saturation,
banners_matrix=banners_carryover_matrix,
original_spends=data,
N=N
)
# Solve optimization
print("Starting optimization with scipy.optimize.minimize (SLSQP)...")
solution = minimize(
fun=partial_objective,
x0=initial_channel_totals,
bounds=bounds,
method="SLSQP",
jac="3-point",
options={
"maxiter": 1000,
"disp": True,
"ftol": 1e-9
},
constraints={
"type": "eq",
"fun": budget_constraint,
"args": (original_total_budget,)
}
)
# Extract optimized channel totals
tv_total_optimized, radio_total_optimized, banners_total_optimized = solution.x
# Display results
print("\n" + "="*60)
print("SCIPY OPTIMIZATION RESULTS (SLSQP)")
print("="*60)
print(f"Status: {solution.message}")
print(f"Success: {solution.success}")
print(f"\nOptimal Total Sales: ${-solution.fun:,.2f}") # Negative because we minimized
print(f"\n--- Budget Allocation ---")
print(f"Original Total Budget: ${original_total_budget:,.2f}")
print(f"Optimized Total Budget: ${solution.x.sum():,.2f}")
print(f"\n--- Spend per Channel ---")
print(f"TV: ${tv_total_optimized:,.2f} (Original: ${data['TV'].sum():,.2f})")
print(f"Radio: ${radio_total_optimized:,.2f} (Original: ${data['Radio'].sum():,.2f})")
print(f"Banners: ${banners_total_optimized:,.2f} (Original: ${data['Banners'].sum():,.2f})")
print(f"\n--- Sales Improvement ---")
print(f"Actual Sales (from data): ${y_actual.sum():,.2f}")
print(f"Optimized Sales (SLSQP): ${-solution.fun:,.2f}")
print(f"Improvement: ${-solution.fun - y_actual.sum():,.2f}")
print(f"Improvement %: {((-solution.fun - y_actual.sum()) / y_actual.sum() * 100):.2f}%")
# Compare with CVXPY results
print(f"\n--- Comparison with CVXPY ---")
print(f"CVXPY Sales: ${problem.value:,.2f}")
print(f"SLSQP Sales: ${-solution.fun:,.2f}")
print(f"Difference: ${abs(problem.value - (-solution.fun)):,.2f}")
Starting optimization with scipy.optimize.minimize (SLSQP)...
Optimization terminated successfully (Exit mode 0)
Current function value: -2237678.352585092
Iterations: 32
Function evaluations: 224
Gradient evaluations: 32
============================================================
SCIPY OPTIMIZATION RESULTS (SLSQP)
============================================================
Status: Optimization terminated successfully
Success: True
Optimal Total Sales: $2,237,678.35
--- Budget Allocation ---
Original Total Budget: $1,336,103.05
Optimized Total Budget: $1,336,103.05
--- Spend per Channel ---
TV: $551,512.71 (Original: $589,241.53)
Radio: $93,956.04 (Original: $442,717.01)
Banners: $690,634.30 (Original: $304,144.51)
--- Sales Improvement ---
Actual Sales (from data): $2,133,628.30
Optimized Sales (SLSQP): $2,237,678.35
Improvement: $104,050.05
Improvement %: 4.88%
--- Comparison with CVXPY ---
CVXPY Sales: $2,482,938.10
SLSQP Sales: $2,237,678.35
Difference: $245,259.74
Comparison of Budget Allocation Methods¶
In [ ]:
# @title
import matplotlib.pyplot as plt
import numpy as np
# Prepare data for comparison
methods = ['Original', 'CVXPY', 'SLSQP']
# Budget allocations
tv_spends = [data['TV'].sum(), tv.value.sum(), tv_total_optimized]
radio_spends = [data['Radio'].sum(), radio.value.sum(), radio_total_optimized]
banners_spends = [data['Banners'].sum(), banners.value.sum(), banners_total_optimized]
# Total sales
sales_values = [y_actual.sum(), problem.value, -solution.fun]
# Create comparison plots
fig, axes = plt.subplots(2, 2, figsize=(15, 10), dpi=150)
# 1. Budget Allocation Comparison (Stacked Bar)
x = np.arange(len(methods))
width = 0.5
axes[0, 0].bar(x, tv_spends, width, label='TV', color='#1f77b4')
axes[0, 0].bar(x, radio_spends, width, bottom=tv_spends, label='Radio', color='#ff7f0e')
axes[0, 0].bar(x, banners_spends, width, bottom=np.array(tv_spends)+np.array(radio_spends),
label='Banners', color='#2ca02c')
axes[0, 0].set_ylabel('Budget ($)', fontsize=11)
axes[0, 0].set_title('Budget Allocation by Method', fontsize=12, fontweight='bold')
axes[0, 0].set_xticks(x)
axes[0, 0].set_xticklabels(methods)
axes[0, 0].legend()
axes[0, 0].grid(axis='y', alpha=0.3)
# 2. Total Sales Comparison (Bar)
colors = ['#666666', '#d62728', '#9467bd']
bars = axes[0, 1].bar(methods, sales_values, color=colors, alpha=0.8, edgecolor='black')
axes[0, 1].set_ylabel('Total Sales ($)', fontsize=11)
axes[0, 1].set_title('Total Sales by Method', fontsize=12, fontweight='bold')
axes[0, 1].grid(axis='y', alpha=0.3)
# Add value labels on bars
for bar, val in zip(bars, sales_values):
height = bar.get_height()
axes[0, 1].text(bar.get_x() + bar.get_width()/2., height,
f'${val:,.0f}',
ha='center', va='bottom', fontsize=9)
# 3. Channel Budget Comparison (Grouped Bar)
x = np.arange(len(methods))
width = 0.25
axes[1, 0].bar(x - width, tv_spends, width, label='TV', color='#1f77b4', alpha=0.8)
axes[1, 0].bar(x, radio_spends, width, label='Radio', color='#ff7f0e', alpha=0.8)
axes[1, 0].bar(x + width, banners_spends, width, label='Banners', color='#2ca02c', alpha=0.8)
axes[1, 0].set_ylabel('Spend ($)', fontsize=11)
axes[1, 0].set_title('Channel Spend Comparison', fontsize=12, fontweight='bold')
axes[1, 0].set_xticks(x)
axes[1, 0].set_xticklabels(methods)
axes[1, 0].legend()
axes[1, 0].grid(axis='y', alpha=0.3)
# 4. Sales Improvement Percentage (Bar)
improvement_cvxpy = ((problem.value - y_actual.sum()) / y_actual.sum()) * 100
improvement_slsqp = ((-solution.fun - y_actual.sum()) / y_actual.sum()) * 100
improvements = [0, improvement_cvxpy, improvement_slsqp]
colors_imp = ['#666666', '#d62728', '#9467bd']
bars = axes[1, 1].bar(methods, improvements, color=colors_imp, alpha=0.8, edgecolor='black')
axes[1, 1].set_ylabel('Sales Improvement (%)', fontsize=11)
axes[1, 1].set_title('Sales Improvement vs Original', fontsize=12, fontweight='bold')
axes[1, 1].axhline(y=0, color='black', linestyle='-', linewidth=0.8)
axes[1, 1].grid(axis='y', alpha=0.3)
# Add value labels on bars
for bar, val in zip(bars, improvements):
height = bar.get_height()
axes[1, 1].text(bar.get_x() + bar.get_width()/2., height,
f'{val:.2f}%',
ha='center', va='bottom' if val > 0 else 'top', fontsize=9)
plt.tight_layout()
plt.show()
# Print detailed comparison table
print("\n" + "="*80)
print("DETAILED COMPARISON: CVXPY vs SLSQP")
print("="*80)
print(f"{'Metric':<30} {'CVXPY':<20} {'SLSQP':<20} {'Difference':<15}")
print("-"*80)
print(f"{'Total Sales':<30} ${problem.value:>15,.2f} ${-solution.fun:>15,.2f} ${abs(problem.value - (-solution.fun)):>12,.2f}")
print(f"{'TV Spend':<30} ${tv.value.sum():>15,.2f} ${tv_total_optimized:>15,.2f} ${abs(tv.value.sum() - tv_total_optimized):>12,.2f}")
print(f"{'Radio Spend':<30} ${radio.value.sum():>15,.2f} ${radio_total_optimized:>15,.2f} ${abs(radio.value.sum() - radio_total_optimized):>12,.2f}")
print(f"{'Banners Spend':<30} ${banners.value.sum():>15,.2f} ${banners_total_optimized:>15,.2f} ${abs(banners.value.sum() - banners_total_optimized):>12,.2f}")
print(f"{'Sales Improvement %':<30} {improvement_cvxpy:>15.2f}% {improvement_slsqp:>15.2f}% {abs(improvement_cvxpy - improvement_slsqp):>12.2f}%")
print("="*80)
================================================================================ DETAILED COMPARISON: CVXPY vs SLSQP ================================================================================ Metric CVXPY SLSQP Difference -------------------------------------------------------------------------------- Total Sales $ 2,482,938.10 $ 2,237,678.35 $ 245,259.74 TV Spend $ 609,141.90 $ 551,512.71 $ 57,629.20 Radio Spend $ 0.08 $ 93,956.04 $ 93,955.96 Banners Spend $ 726,961.06 $ 690,634.30 $ 36,326.76 Sales Improvement % 16.37% 4.88% 11.49% ================================================================================
CVXPY Constraints:¶
Total: 601 constraints Budget constraint type: Inequality (≤) - can spend up to the budget Variable count: 600 variables (200 time periods × 3 channels)
scipy.optimize.minimize Constraints:¶
Total: 3 bound constraints + 1 equality constraint = 4 constraints Budget constraint type: Equality (=) - must spend exactly the budget Variable count: 3 variables (total per channel)